function [I] = woodhull(v, Vr, lam)
      % Woodhull fit from Washburn et al, J. Neurosci, 1998. 
      % input is IV for current
      % parameters are:
      % lam(1) is g scaling factor
      % lam(2) is blocker concentration
      % lam(3) is Kd(0)
      % lam(4) is z*(1-delta)
      % alpha is RT/F (must be fixed when called!)
      % Vr is reversal potential...
      % G0 is set to be the function given in that paper...
      alpha = 25.3;
      g0 = 0.3 + 0.6*exp((v-50)/40);
      gv = g0./(1 + (lam(2)./(lam(3)*exp(-lam(4)*v/alpha))));

      I = lam(1).*(v-Vr).*gv;
        plot(v, gv);
        return;
        